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ABSTRACT: The standard picture of electroweak baryogenesis requires slowly expanding 
bubbles. This can be difficult to achieve if the vacuum expectation value (VEV) of a 
gauge singlet scalar field changes appreciably during the electroweak phase transition. It 
is important to determine the bubble wall velocity in this case, since the predicted baryon 
asymmetry can depend sensitively on its value. Here, this calculation is discussed and il¬ 
lustrated in the real singlet extension of the Standard Model. The friction on the bubble 
wall is computed using a kinetic theory approach and including hydrodynamic effects. Wall 
velocities are found to be rather large (v w > 0.2) but compatible with electroweak baryo¬ 
genesis in some portions of the parameter space. If the phase transition is strong enough, 
however, a subsonic solution may not exist, precluding non-local electroweak baryogenesis 
altogether. The results presented here can be used in calculating the baryon asymmetry 
in various singlet-driven scenarios, as well as other features related to cosmological phase 
transitions in the early Universe, such as the resulting spectrum of gravitational radiation. 
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1 Introduction 


Cosmological phase transitions in the early Universe are interesting for a variety of rea¬ 
sons. They can produce observable gravitational radiation [1-4], seed primordial magnetic 
fields [5, 6], affect the abundance of thermal relics [7, 8], and otherwise play an important 
role in the cosmological history of the Universe [9]. Perhaps most notably, such phase 
transitions can give rise to a viable mechanism of baryogenesis, provided the transition is 
first-order. 

A first-order phase transition can occur in the early Universe when two vacua of the 
theory coexist for some range of temperatures. If this is the case, an energy barrier exists 
between the two and the system can transition to the state with lower free energy via 
quantum tunneling or thermal fluctuations [10-13]. Physically, this corresponds to the 
formation of spherical bubbles in the ambient metastable vacuum. These bubbles grow and 
can reach a steady state expansion velocity. Inside the bubble, some subset of the scalar 
fields have non-zero vacuum expectation values (VEVs), while outside they do not. If the 
non-vanishing condensate breaks the SU(2)l gauge symmetry of the Standard Model (SM), 
as at the electroweak phase transition (EWPT), non-perturbative sphaleron transitions, 
which violate B + L , will be quenched inside the bubble and active outside. These processes 
can act on chiral charge currents diffusing in front of the wall to source a net baryon 
asymmetry. If the sphaleron rate is significantly suppressed in the broken electroweak 
phase, the asymmetry can be frozen in once captured by the expanding bubble. Roughly, 
this requires [14, 15, 116] 


(h) 

T 


> 1 , 


( 1 . 1 ) 


where T is a temperature associated with the phase transition and h is some combination 
of fields charged under SU (2 )l usually identified with the Standard Model-like Higgs. This 
picture is known as "non-local" or "transport-driven" electroweak baryogenesis (EWB) and 
is an elegant explanation for the origin of the observed baryon asymmetry of the Universe 
[17-22], 

Clearly this mechanism relies on several requirements beyond sphaleron suppression 
in the broken phase. One must ensure a significant amount of CP-violation to source the 
chiral charge currents in front of the wall. This is well-known, has been studied exten¬ 
sively, and has motivated several experiments in search of CP-violating signatures, such as 
permanent electric dipole moments [23]. There is another, somewhat less-appreciated, re¬ 
quirement for successful baryogenesis through this mechanism: bubbles must expand slowly 
enough for sphalerons to convert a significant fraction of the CP-asymmetry to a net baryon 
density [19-21, 24], In other words, the diffusion of the chiral plasma excitations must be 
efficient in front of the bubble wall 1 . 

Precisely how slowly the wall must move in this scenario depends on several factors, 
including the amount of CP-violation and the details of diffusion in front of the bubble. 


1 Bubbles must also not expand too slowly; otherwise a quasi-equilibrium situation is reached and the 
net baryon density is equilibrated away. This is typically not a problem in baryogenesis scenarios, since it 
requires bubbles moving very slowly, with v w < 0.01 [19] 
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However, it is generally the case that bubble expansion must at least be slower than the 
speed of sound in the plasma 2 , c s ~ l/\/3 ~ 0.58. Otherwise, diffusion in front of the 
wall will be very inefficient. Even if the wall moves subsonically, the predicted value of 
the steady-state wall velocity, v w , is an important input into any non-local electroweak 
baryogenesis calculation. Previous studies have found that the predicted baryon asymmetry 
typically peaks around v w ~ 0.01, falling off as ~ l/v w or faster for larger values [26-29] (for 
velocities much smaller than 0.01, the sphalerons can begin to equilibrate the asymmetry). 
In some cases the dependence of the predicted asymmetry on the wall velocity can be less 
severe [30-32], as the scaling hinges on the form of the dominant CP-violating source (see 
e.g. Refs. [33-35]). Nevertheless, a determination of the wall velocity is often important 
even for rough estimates of the baryon asymmetry when the wall is expected to move 
quickly. 

The electroweak phase transition in the Standard Model is not first-order [36, 37]. Thus, 
electroweak baryogenesis necessarily requires some physics beyond the SM. There have been 
many such scenarios proposed in the literature. One of the most popular and straightfor¬ 
ward involves augmenting the Standard Model Higgs sector by a gauge singlet scalar field 
[38-51], providing a tree-level cubic term in the effective potential. Such a cubic term can 
easily give rise to the barrier required for a first-order transition, and the singlet nature 
of the new state(s) in many cases can ensure compatibility with current phenomenological 
constraints [47-51], including the observation of the 125 GeV Standard Model-like Higgs at 
the LHC [52, 53]. This class of models, though popular and simple, is often expected to 
produce fast-moving bubble walls when the singlet field VEV changes appreciably during 
the electroweak phase transition [44], This is simply because the additional field direction 
contributes to the pressure difference between the phases, which drives the expansion of the 
bubble, but does not experience a substantial drag force from the plasma 3 . Although this 
fact was recognized several years ago [44], studies of EWB typically focus on the strength of 
the phase transition and assume C-P-violation can be added in separately without consid¬ 
ering the effects of the bubble wall dynamics on generating the baryon asymmetry. In fact, 
while many studies have since considered baryogenesis in these scenarios with a changing 
singlet VEV, the wall velocity has never been directly calculated 4 , nor has it been shown 
that the resulting bubble walls can propagate subsonically as required for successful EWB. 
Our aim here is to fill this gap. 

The bubble wall velocity is an important quantity to compute apart from baryogenesis 
considerations. For example, models with additional gauge singlets that predict a strong 

2 More precisely, it is the wall velocity relative to the fluid in front of the bubble that should be sub¬ 
sonic [25]. However, in what follows the fluid velocity in the symmetric phase will always be perturbatively 
small, and so we will simply require the wall velocity to be subsonic in the rest frame of the fluid far from 
the bubble. For further discussion on this point, see Ref. [25]. 

2 If the singlet field is instead approximately stabilized during the electroweak phase transition, bubbles 
can expand significantly more slowly [54]). 

4 In Ref. [48], we performed a very rough estimate of the wall velocity in the NMSSM. Several important 
terms in the Boltzmann equations were dropped, likely resulting in a significant under-estimate of the wall 
velocity. A full microphysical calculation of v w in the NMSSM does not currently exist in the literature. 
Such a study can be undertaken with the methods discussed here and is currently in progress. 
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first-order transition can source gravitational waves through bubble collisions and turbu¬ 
lence (see e.g. Refs. [1, 4, 55-68]). These scenarios may be effectively probed by upcoming 
gravitational wave experiments, such as eLISA [69], or Big Bang Observer [70]. However, 
in order to connect these observations to an underlying theory, one must be able to reliably 
calculate the velocity of the expanding bubble, as well as the other properties of the phase 
transition. 

A detailed calculation of the wall velocity is rather involved. Building on previous 
work [26, 71-75], Moore and Prokopec were the first to calculate the velocity for the Stan¬ 
dard Model case microphysically in Refs. [76, 77] (Ref. [54] also recently revisited this cal¬ 
culation). Five years later, John and Schmidt performed an analogous study in the minimal 
supersymmetric Standard Model (MSSM) with a light scalar top quark (stop) [78]. Around 
the same time, the effects of infrared gauge bosons on the wall velocity were calculated in 
Ref. [79]. Most recently, Ref. [54] extended the results of Moore and Prokopec to other 
SM-like scenarios, including that in which the VEV of a singlet scalar is approximately 
stabilized during the EWPT. To date, these remain the only full microphysical calculations 
of the wall velocity existing in the literature 5 6 . Recent years have seen progress in match¬ 
ing models onto these existing results [54, 80-82] and in the hydrodynamic considerations 
associated with bubble wall expansion (see e.g. [83-89]). 

Scenarios in which an the VEV of an additional singlet scalar field changes appreciably 
during the transition merit separate consideration 5 . This is because one must account for 
the friction on the singlet field. Neglecting these contributions can lead to a drastic over¬ 
estimate of the wall velocity. A proper treatment requires computing several new classes 
of interaction rates in the plasma, which can be rather involved. Also, the additional 
field direction complicates the equations of motion for the condensates. Nevertheless, the 
calculation can be done and is especially important given the current status of electroweak 
baryogenesis in light of collider searches. Until recently, the MSSM light stop scenario 
[90, 91] was considered by many to be the most plausible setting for electroweak baryogenesis 
beyond the SM. Now light stops are in severe tension with both direct LHC searches [92, 93] 
and measurements of the Standard Model-like Higgs couplings [94-97]. Similar conclusions 
hold true for many different models relying on large thermally-induced cubic terms to 
strengthen the phase transition [94, 98, 99]. This situation has led to a renewed interest in 
singlet-driven scenarios, since they can be much more elusive at colliders [50]. An analysis 
of the wall dynamics would mark an important step forward in understanding electroweak 
baryogenesis in these models. 

5 By ‘microphysical’, we mean calculations explicitly computing the friction exerted by the plasma on the 
wall, as opposed to those using a phenomenological viscosity parameter. The former involves determining the 
various interaction rates in the electroweak plasma and solving for the deviations from thermal equilibrium 
around the bubble wall, as we describe in detail below. 

6 In the remainder of this study we will take ‘singlet-driven scenarios’ to refer to those in which the 
singlet VEV changes non-negligibly during the transition. This can occur even in models with a discrete Z 2 
symmetry at T = 0. Examples in which the singlet VEV is approximately stabilized during the EWPT can 
be treated by the techniques developed for the Standard Model (or MSSM, if the singlet contributions to the 
finite-T cubic term are large) since only the Higgs field is involved in the transition (see e.g. Refs. [54, 82] 
for an application of this approach). 
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The goal of this study will be to demonstrate how the electroweak bubble wall velocity 
can be calculated and to extract some general features of the wall dynamics in singlet-driven 
scenarios . We will generally follow the strategy and techniques developed in Refs. [54, 76- 
79] but modified to account for the singlet field direction. The methods described will be 
applicable to many different models, although for the sake of simplicity we will frame our 
discussion in the real singlet extension of the SM, sometimes known as the ‘xSM’ [43]. This 
scenario should encapsulate the most relevant features of models with singlet-driven first- 
order phase transitions. Notably, the xSM does not feature any new sources of CP-violation, 
which could in principle significantly alter the wall dynamics when included. We comment 
further on this below and the reader should keep this in mind as we proceed. 

We will focus on two different schemes for calculating the wall velocity. The first is 
explicitly gauge-independent and neglects the contributions to the effective potential and 
friction from the SU(2)l gauge bosons, while the second includes them. For slower bubble 
walls (such that the friction on the singlet field is large), both calculations yield similar 
results, while for faster walls the gauge boson contributions become increasingly important 
in slowing down the expansion. 

The remainder of this study is structured as follows. In Section 2, we introduce the 
singlet-driven scenario and the finite temperature effective potential which will be used 
throughout our study. We then move on to computing the wall velocity. Following Ref. [77] , 
the calculation can be broken down into several parts. First, the phase transition properties 
must be computed and the temperature near the bubble inferred from hydrodynamic con¬ 
siderations, as discussed in Sec. 3. There we also discuss the equations of motion (EOMs) 
for the bubble wall and consider the simple case of wall velocities in the ultra-relativistic 
limit. Possible values for the steady state wall velocity are those such that the equations 
of motion for the wall are satisfied. The EOMs depend on the deviations from thermal 
equilibrium of all plasma excitations in front of the wall. These are discussed in Sec. 4. We 
next move on to solving the system of equations for the deviations from equilibrium and the 
equations of motion in Sec. 5. The calculation is then applied to the parameter space of the 
xSM consistent with all phenomenological constraints in Sec. 6. We find that sufficiently 
strong phase transitions may possess no subsonic solutions, and that v w > 0.2 for points 
with (h)/T n > 1 in the parameter regions considered. We infer that bubbles may expand 
slowly enough for singlet-driven electroweak baryogenesis, but only in certain portions of 
the parameter space. Results for the bubble wall profiles are also presented. Finally, our 
main findings and conclusions are summarized in Sec. 7. We also provide a brief appendix 
which compares the interaction rates we have calculated with those appearing elsewhere in 
the literature. 

Before proceeding it is important to note that there are other electroweak baryogenesis 
scenarios that do not rely on diffusion in front of the wall, and hence that do not require 
slow bubble walls. Local EWB [100-103], in which the baryon number and CP -violation 
occur in the same region at the bubble wall boundary, is one such example'. However, this 
typically leads to a highly suppressed total baryon asymmetry relative to the non-local case, 

7 Cold Electroweak Baryogenesis [104-106] also falls into this category. 
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since the sphaleron transitions turn off near the outer edge of the bubble wall [19]. More 
recently, an interesting scenario was presented in Ref. [107], in which bubbles can expand 
quickly enough to significantly reheat the plasma inside the bubble. Secondary bubbles can 
then nucleate near which transport-driven baryogenesis can occur. This is an intriguing 
possibility, however it requires a substantial amount of reheating which only occurs for very 
strong phase transitions. The reader should bear these alternative scenarios in mind as we 
proceed. 


2 A Singlet-Extended Higgs Sector 


To study the dynamics of singlet-driven electroweak phase transitions, we will work in 
the real singlet extension of the Standard Model. This simple scenario has been studied 
in depth in the literature, from the standpoint of electroweak baryogenesis, dark matter, 
LHC signatures, and more (see e.g. Refs. [42, 43, 45, 49, 108-113] and references therein). 
The arguments and methods we discuss here can be straightforwardly applied to other 
singlet extensions of the SM, such as the next-to-minimal supersymmetric Standard Model 
(nMSSM [31, 32, 41] or NMSSM [32, 38-40, 47, 48]) and other scenarios with real or complex 
gauge singlets [51]. 

The tree-level potential is taken to be 8 


Vo(H, S)=- ^ 2 {H^H) + X(H^H) 2 + ^ ai (H^H)S + ^a 2 {H^H)S 2 
+ S 2 + -63 S ' 3 + - 64 .S 4 


( 2 . 1 ) 


where S is a real scalar singlet under the Standard Model gauge groups and H is a complex 
SU(2)l doublet. Both the singlet and CP-even neutral component of H are assumed to 
obtain vacuum expectation values during electroweak symmetry breaking. Throughout our 
discussion, we will also assume that both VEVs vanish in the high-temperature phase for 
simplicity (this is discussed further below). Non-vanishing VEVs correspond to minima of 
the effective potential for non-zero background field values 4>h, (j) s . These classical fields are 
those relevant for computing the properties of the phase transition. At a given temperature, 
we can expand H and S about the background fields, 


H t 


(t>h{T) + h + ic/) 0 \ 

V2 )' 


S = MT) + s. 


( 2 . 2 ) 


At zero temperature, (j>h(T = 0) = v = 246 GeV. The zero temperature singlet VEV, 
4> S (T = 0 ) = v s can vary. 

Throughout our study we will identify h with the Standard Model-like Higgs discovered 
at the LHC [52, 53] and take s to be a pure singlet with no mixing at tree-level. The 
phenomenology of this setup and our choices for the various parameters are detailed in 
Section 6.1 below. 

8 One is free to shift the singlet field value such that the T = 0 tadpole is removed [45]. 
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2.1 The Effective Potential 


For a homogeneous background field configuration (j)(x) = 4>, the ground state of the theory 
corresponds to a minimum of the effective potential V e f At one loop, 14ff is given by 
the tree-level potential (expanded around the background fields), modified by additional 
Coleman-Weinberg terms. 

At finite temperature and density the physical ground state of the theory is altered by 
the interactions of the scalar field 0 with the ambient plasma. The vacua of the theory can 
then be determined from the finite-temperature effective potential, V e s(4>, T ). In the simple 
case involving one background held, it is given by [14] 


V t {^T) 


rp 4 
27T 2 


Y, ±NiJ ± 



(2.3) 


where the plus and minus signs correspond to the bosonic and fermionic contributions, 
respectively, and the N t are the associated number of degrees of freedom for the species 
i. This expression generalizes straightforwardly to the case of more than one background 
held. The functions J± are given by 


J±(x) = / dyy 2 log 1 =F exp(-^x 2 + y 2 ) 

Jo 1 


(2.4) 


In the high-temperature limit they admit a useful expansion, given by 


4 f m 2 \ 7t 4 T 4 ir 2 m 2 T 2 Tir(m 2 ) 3 / 2 (m 4 ) m 2 

+ ^^2 J = “ 45 + 12 6 32“ ° g a b T 2 ’ 

4 (m 2 \ _ 7 tt 4 T 4 n 2 m 2 T 2 (m 4 ) m 2 

- J 360 24 32~ ° g a/T 2 ’ 


with a b = 16 7 r 2 e 3//2 ~ 27s , a/ = 7 r 2 e 3 / 2 ^ 27B , and 7 ^ the Euler-Mascheroni constant. Note 
that the thermal contributions above correspond to momentum integrals of equilibrium 
distribution functions for all species in the plasma coupled to (f> [14, 44], 


2.2 Gauge-Invariance 

The finite temperature effective potential is only gauge-invariant at its extrema [114, 115]. 
Thus, tunneling calculations depending on the potential away from the local minima are 
in general gauge-dependent. This will result in a gauge-dependent determination of the 
nucleation temperature for the phase transition, T n , and ultimately the wall velocity. To 
avoid this as much as possible, our primary analysis will only consider terms in the effec¬ 
tive potential which are explicitly gauge-invariant. Thus, we will not include the T = 0 
Coleman-Weinberg corrections, or the finite temperature cubic and tadpole terms in the 
high-temperature effective potential (gauge-dependence in the tadpole may enter at higher 
perturbative order [49]). This is precisely the strategy followed by Ref. [49] in analyzing the 
phase transition properties of the xSM. The finite-temperature effective potential in this 
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case becomes 


Vesi^hi <t>siT) — — 2^ 2 4 >\ + + ^ a l (frs&h + 


b\T ' 2 


ti? 2 


H—(9^2 + 3gf + 12 y^ + 24A + 202 ) H—(2a2 + 364 ). 


24 


( 2 . 6 ) 

Although morally satisfying, the gauge-invariant approach has the disadvantage of 
sometimes neglecting numerically important contributions to the effective potential. While 
it should capture the physics we are interested in, namely the singlet contributions to the 
potential and the friction on the expanding bubble, it is important to consider the effects of 
the neglected terms, especially the gauge boson cubic term 9 . To this end, we will also show 
results for the wall velocities including the gauge boson cubic term and friction in Landau 
gauge. In this case the effective potential in Eq. 2.6 is modified by the additional term 

k lT 


a t /-cubic 




12 vr 


4 92 + g {d 2 + 9i) 


3/2 


(2.7) 


We will find that these contributions (and the friction from gauge bosons) are numerically 
significant in some cases, especially for faster moving bubble walls. 


3 Preliminaries: Phase Transitions, Hydrodynamics, and the Wall Equa¬ 
tions of Motion 

For a given point in the model parameter space with a first-order electroweak phase transi¬ 
tion, we are interested in determining the steady-state velocity of the bubble wall separating 
the electroweak-symmetric and broken phases. To do so we must first determine the prop¬ 
erties of the phase transition (most importantly its characteristic temperature), as well as 
the equations of motion for the scalar field condensates. Our discussion will roughly follow 
that of Ref. [77], which we draw from frequently throughout the remainder of this study. 

3.1 Bubble Nucleation 

Using the effective potential in Eq. 2.6, first-order transitions can occur when two local 
minima coexist for some range of temperatures. The background fields can then “tunnel" 
from the origin to the new vacuum, in which (j)^, <j) s ^ 0. This can begin to occur below 
the critical temperature, T c , at which the two relevant vacua are degenerate. Bubbles 
begin to nucleate efficiently at the nuceation temperature, T n . determined by requiring the 
expectation value for one bubble to nucleate per Hubble volume to be ~ 0(1). At finite 
temperature, the nucleation probability is determined by the 0(3)-symmetric instanton 
interpolating between the metastable and and true vacua with the lowest ratio of three- 
dimensional Euclidean action 10 to temperature, S 3 /T [12-14], The tunneling probability 

9 The inclusion of the tadpole term acts primarily to shift the high-temperature minimum away from 
( j>„ = 0. We have computed the wall velocities for several scenarios with this term included and obtain 
values similar to those found neglecting the tadpole. For simplicity we will not consider this term further 
in this study, although our methods can be straightforwardly modified to include it. 

10 For all the cases we consider, the nucleation temperatures are much larger than the inverse radii of the 
instantons, and so the 0(3) bounce is indeed the relevant quantity to consider. 
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per unit volume is then given by 


y = A(T)e~ s ^ T . (3.1) 

where the pre-factor A(T ) is only weakly temperature-dependent. Using dimensional anal¬ 
ysis to estimate A(T ) and assuming typical electroweak temperatures, one hnds that the 
nucleation temperature is approximately determined by Ss(T n )/T n ~ 140 [14], We adopt 
this definition in the rest of our study. 

Inside the nucleated bubble, both the VEV of the Higgs and singlet fields will be non¬ 
zero by assumption. From the standpoint of electroweak baryogenesis, only the value of tph 
is important for sphaleron suppression. We will therefore define a strongly first-order phase 
transition, occurring at the nucleation temperature T n , by [116] 

> 1 - (3 

This condition is free from explicit gauge-dependence in our primary setup, since we have 
neglected all gauge-dependent terms in the effective potential. It should be noted that 
the baryon number preservation condition above still contains several implicit assumptions, 
discussed in detail in Ref. [15]. 

3.2 Temperature Variations 

The nucleation temperature defined above is that of the ambient plasma when bubbles 
begin to form efficiently. Once formed, however, the temperature is no longer homogeneous. 
The phase transition releases latent heat into the plasma and the expansion of a subsonic 
bubble heats up the medium in front of it. The temperature in the broken phase will 
thus differ from that immediately outside the bubble, which in turn is not the same as the 
typical nucleation temperature of the bubble. To relate these various quantities requires 
a treatment of the plasma hydrodynamics. These changes in temperature can have large 
effects on the expansion of the bubble [77], and so we must take them into account. 

Far away from the bubble, the relevant temperature is that at which bubble nucleation 
occurs, T n . We wish to obtain the temperature in the vicinity of bubble wall. To do so, 
let us consider the wall-plasma system, with the plasma modeled as a perfect relativistic 
fluid. Hydrodynamic equations can be obtained by requiring conservation of the wall-fluid 
stress-energy tensor [71], 



a,T>‘- = = o. (3.3) 

We define the ‘fluid or ‘plasma frame’ such that the fluid is at rest far from the bubble 
and in its center. This is the frame which we use to define the wall velocity v w and the 
wall profile parameters. Solutions to the fluid equations in the plasma frame can typically 
be classified as either ‘detonations’, in which the bubble velocity exceeds the sound speed 
c s in the plasma, or ‘deflagrations’, in which v w < c s n . Successful subsonic electroweak 

11 There are also ‘hybrid’ cases; see Ref. [83]. 
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baryogenesis typically requires a deflagration solution, since otherwise diffusion in front of 
the bubble is inefficient. We will restrict ourselves to this case. 

Consider an expanding bubble with free energy V e s((f)-, T_) inside and I4ff(^+) T + ) 
immediately outside (‘±’ subscripts will correspond to quantities outside/inside the bubble). 
The equations of state (EoS) for the two phases can be written as 

P± = 3 a ± (T)Ti - e±(T), p ± = a±(T)T:f + e±(T) (3.4) 

where p± and p± are the pressure and energy density of the fluid in either phase, and 

a±(T) = dVeS [ ^ T),T \ e±(T) = U e ff [<i>±(T),T\ + \a ± (T)T\ (3.5) 

The above form for the equations, taken from Ref. [83], are inspired by the so-called ‘Bag 
EoS’, but involves the temperature-dependent quantities a±(T), e±(T). Fortunately, we can 
safely neglect the temperature dependence in a±. e±. using their values at T = T n . This is 
because the free energy (and hence a±. e±) are dominated by light degrees of freedom, which 
contribute a constant term to the free energy in each phase that does not vary significantly 
between T n and T c for the cases we consider. We find that using a±(T n ), e±(T n ) in Eq. 3.5 
reproduces the full result for the pressure and energy density within a few percent. This is 
fortunate as it allows us to avoid several issues arising for more complicated temperature 
dependence in the EoS, such as the variation of the sound speed in the plasma [87]. 

Ultimately we will take T eq (x ) = T+-\-5Tb g (x) to be the (space-time-dependent) temper¬ 
ature entering the equilibrium distribution functions for the various particles in the plasma. 
We thus need to determine T + , the temperature just outside the bubble. The pressure and 
energy density can be related to the fluid velocities on either side of the phase boundary 
by integrating Eq. 3.3 across the wall. This yields expressions for the fluid velocities v± in 
the wall frame which depend on T±: 


V + V- 


P+-P- 

P+-P- 


v+ = p- + p + 

p + + p- ' 


(3.6) 


These velocities can be simply transformed to their analogs in the fluid frame v± via v± = 
v w — \v±\/(l — u^jn-i-l). Note that for a subsonic deflagration, V- = 0 and so v w = — u_, as 
discussed in e.g. Refs. [71, 83]. 

One then needs to relate the temperature T + to T n . In the subsonic deflagration case, 
the bubble wall is preceded by a shock front moving with velocity v s h in the fluid frame. 
The temperatures T \ i 2 and fluid velocities v\p, on either side of the shock front (in its rest 
frame) will be different; the equation of state is however the same (again neglecting small 
temperature variations in a+(T), e+(T)). We will use the subscripts 1,2 to denote quantities 
inside and outside the shock front, respectively. One can again integrate across the interface 
and use the fact that the fluid is at rest beyond the shock front (i.e. V 2 = 0 —> V 2 = —v s h) 
with temperature T 2 = T n . This yields an expression for v\ in terms of T\ and T n : 


2 _ 3 n+Tt 

Vl 9 Tf + 3T4' 


(3.7) 
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Figure 1 . Illustration of the competing forces acting on the bubble wall that ultimately determine 
v w . The steady state wall velocity is such that the vacuum energy difference between the phases 
(AVt=o) is balanced by the friction provided by the interactions of the wall with the plasma. 


Again, the corresponding fluid velocity in the fluid frame v\ is simply given by velocity 
addition, v\ = 1 . 

Throughout our calculation we neglect the curvature of the bubble wall. In this ap¬ 
proximation, the temperatures and fluid velocities (in the fluid frame) between the bubble 
wall and the shock wave are simply constant [84], and so one can set 


v\{Ti,T n ) « v + (T + , v w ), Ti « T + 


(3.8) 


and solve for T + in terms of T n , v w . Previous studies suggest that using the planar approxi¬ 
mation instead of the full solutions to the spherical hydrodynamic equations can reproduce 
the full result for the wall velocity to within a few percent [82]. 

With the temperature T + and the static properties of the phase transition determined 
in this way, we can now consider the asymptotic behavior of the bubble after its formation. 


3.3 Wall Equations of Motion 

The main object for our analysis will be the bubble wall equations of motion correspond¬ 
ing to the set of scalar fields fa = < j>h,(t>s- These can be derived by requiring conservation 
of the energy-momentum tensor for the scalar held condensates computed in a WKB ap¬ 
proximation [77], or directly from the Kadanoff-Baym equations [54], We are interested in 
the stationary limit of the equations of motion in the plasma frame; that is, we want to 
investigate the bubble wall once it has reached its terminal velocity (if it exists), with the 
pressure driving the expansion precisely counterbalanced by the drag force exerted on the 
bubble by the plasma. This is illustrated in Fig. 1. 

Neglecting the curvature of the bubble, in the rest frame of a stationary (non-accelerating) 
bubble wall all functions will be depend only on z. the distance from the phase boundary. 
Consequently, in the plasma frame, all functions depend only on the coordinate x = z+v w t, 
where v w is the wall velocity in the plasma frame and we have assumed that the wall is 
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moving to the left. In the stationary wall limit, the equations of motion then simplify to 


-(! - v l 


b „ + dV{<lH,T) 


d<t>i 


x - dm 2 -(fa 
j 


d 3 p 

(27t) 3 2£1 


-5fj(p,x) = 0 


(3.9) 


where primes indicate differentiation with respect to x. Here the sum is over all fields 
coupling to the scalar held fa. Ej is the (space-time-dependent) energy of the particle j. 
Ej = yjp 2 + m 2 (x), and 6 fj is the deviation from the equilibrium distribution function for 
the species j. 

Solutions to the above equations of motion typically only exist for one subsonic value 
of the constant v w . This is the quantity we wish to determine. To do so, one must find 
profiles fa(x) such that Eq. 3.9 is satisfied, which in turn requires solving for the deviations 
from equilibrium of the various species in the plasma. These deviations, along with the 
equilibrium contributions, are responsible for the drag force on the bubble wall. Unfortu¬ 
nately, the 5fj depend non-trivially on v w and the bubble profile, so Eq. 3.9 represents a 
set of integro-differential equations. 


3.4 Aside: Runaway Bubbles and Tree-level Cubic Terms 

Before moving on to the case of non-relativistic bubbles (relevant for electroweak baryo- 
genesis), we can begin by considering the wall dynamics in a simple limit: that of ultra- 
relativistic, “runaway” bubbles [44], with Lorentz factor 7 > 1. In this case, the friction on 
the bubble from the plasma in the large -7 limit is too small to counterbalance the pressure 
difference between the vacua, which drives the expansion. Ref. [44] showed that this situ¬ 
ation is common in singlet-driven transitions, so it is important to review this case before 
moving on to the non-relativistic regime. 

Following Ref. [44], a runaway solution to the equations of motion exists provided 

V eS {T = 0,c/) + ) -V eS (T = 0,<j)_) + [m 2 ((/) + )-m 2 ^)) J ^U^/o,i(p, <P+) > 0 

(3.10) 

at the nucleation temperature. Here, /o is the equilibrium distribution function of the 
species i, and <j>± are the held values at the minima of the potential. In the high-T limit, 
there is a simple interpretation of this criterion in terms of the high-temperature expansion 
of the thermal effective potential: a runaway solution will exist if it is energetically favorable 
to tunnel to the broken phase in the ‘mean-held’ potential, obtained by retaining only the 
T 2 terms in Eq. 2.6. In other words, 

V e g° cublc (<f> + , T n ) > Vfa° cublc ((/>_, T n ) => runaway solution exists. (3-11) 

The above expression indicates that all points found with a first-order phase transition 
in our gauge-invariant approach (retaining only the quadratic hnite-T terms) would fea¬ 
ture an ultra-relativistic wall solution if there were no other contributions to the effective 
potential. This may appear incompatible with our goal of determining subsonic solutions 
to the equations of motion but it is not. First of all, including the finite temperature 
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cubic term inevitably changes the the transition temperature and the effective potential 
at that temperature. This can cause the same parameter space point to instead feature 
Vg cublc (</>+, T n ) < Vg cublc (0_, T n ) , and hence no runaway solution. We indeed find this 
to be the case for most points considered when including the gauge boson cubic term in 
our parameter scans. Even if a runaway solution exists for the EOMs including the full 
finite-T effective potential, there is another important caveat. The criterion in Eq. 3.11 
assumes that the bubble is in the ultra-relativistic regime to begin with. However it is 
instead possible for the friction to prevent the bubble from ever reaching such large veloci¬ 
ties required for Eq. 3.10 to be valid. In fact, hydrodynamic effects alone can obstruct the 
wall from expanding ultra-relativistically [117]. Thus, even if a particular parameter space 
point admits a runaway solution, it may not be realized if a subsonic stationary solution 
exists. On the other hand, even if no runaway solution exists, one with v w > c s might. 
The reader should thus bear in mind that our approach will find subsonic solutions to the 
equations of motion, not guarantee that they are realized. This is also true of previous 
studies [54, 76-78]. 

Equation 3.10 shows that the friction force acting on the wall takes a very simple 
form in the 7 > 1 limit. This is not the case for the subsonic walls we are interested 
in. Determining the wall velocity in the 7 ~ 1 regime requires a careful calculation of the 
various deviations from equilibrium in the plasma. This is what we discuss in the following 
section. 

4 Kinetic Theory and Deviations from Equilibrium 
4.1 Setup 

With the temperature T = T + inferred from hydronamic considerations, the first step 
towards solving the bubble wall equations of motion in the non-relativistic (7 ~ 1 ) case is 
determining the distribution functions for the various excitations appearing in Eq. 3.9. 
To do so, we will primarily utilize a perturbative effective kinetic theory approach [118, 
119], as in previous studies [76-78] (we will take a somewhat different approach for the 
corresponding gauge boson friction, which should be modeled classically as discussed below). 
This treatment applies to weakly coupled excitations with local interactions and short 
wavelengths compared to the length scale of the bubble wall in the plasma frame, i.e. 

E » J— (4-1) 

-L-'w 

where L w is the wall width. Typical momenta are of order p ~ T, but softer excitations 
will be present in the plasma as well. We will assume that the kinetic theory description 
is viable in the range p > gT, which is reasonable for the particles we will be interested in 
given the values we find for the wall widths. Here and throughout this section g represents 
a generic dimensionless coupling of the theory 12 that is assumed to be small. Infrared 

12 The coupling g should be thought of as some combination of couplings entering the thermal and zero 
temperature masses of the particle in question. In other words, we assume parametrically that m ~ gT ~ gif) 
near the electroweak phase transition. 
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(IR) excitations with momenta p<T will not be captured by this treatment, since their 
interactions cannot be properly described by a local collision term. These contributions 
can be important for the bosonic species [79], but the perturbative effective kinetic theory 
should provide an adequate estimate of the damping force on the bubble wall, provided 
that very infrared excitations are equilibrated quickly [77], as we will assume for most of 
the species we are interested in 13 . 

In the effective kinetic theory we consider, the quasiparticle distribution function for 
the species i satisfies the Boltzmann equation 

(4 - 2) 


in the fluid frame, where C[f]i is a local collision integral. The collision term involves all 
interactions of the species i with all other excitations in the plasma. It can be written as 

[119] 




2 Er, 


jmn 


d 3 kd 3 p'd 3 k' 

(2Tr)^2Ei i 2Epi2Eki 


I Mij^mnip, k-,p ', k') | 2 (27 t) 4 6(p + k-p' - k') 


X Vij^mn fj(k ), fmip'), fn{k')\ 

(4.3) 

where the sum is over all 4-body processes ij —> mn, with the momenta labeled as p, p ', A/, 
and k moving clockwise around the diagram starting with particle i. The matrix elements 
include finite-temperature effects (discussed below) and are summed over helicities and 
colors of all four external quasiparticles, then divided by the number of degrees of freedom 
corresponding to species i, Ni ( Nh = 1, Nt = IVj = 6) 14 . The population factor is 


Vij^mn = fifj( 1 ± /m)(l ± fn) ~ fmfn( 1 ± fi)( 1 ± fj) (4.4) 

with the upper (lower) signs corresponding to bosons (fermions) and f a the appropriate 
Bose-Einstein or Fermi-Dirac distribution function for particle a, which we assume to take 
the form 

/„ = (e (B+i * )/T ± l) . (4.5) 

In Eq. 4.3, the prefactor of 1/2 takes care of both the symmetry factor when identical parti¬ 
cles are present in the final state, and the double counting that occurs from the unrestricted 
sum over m and n. 

The Boltzmann equations above apply to all quasiparticles in the plasma satisfying 
Eq. 4.1 with sufficiently high momentum. However, examining Eq. 3.9, we see that only 
the distribution functions of held excitations with significant couplings to the relevant scalar 

lli This was shown to be a poor assumption for the SU(2)l gauge bosons in Ref. [79], which we discuss 
further in Sec. 4.3. Infrared contributions from the Higgs and singlet fields may be important. However, their 
equations of motion are not over-damped as they are for the gauge bosons [79], and so their distributions 
should equilibrate more quickly than those for the gauge fields. 

14 We will neglect any possible CP-violation coupling to the top quark and hence assume that the top 
and anti-top densities are identical. This means we can compute the top perturbations and simply count 
their contribution to the condensate equations of motion twice. 
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fields involved in the phase transition are required. Since these particles have significant 
couplings to the Higgs and singlet scalar fields, we will refer to them as ‘heavy’. Also, 
5fi = 5fi(p,x ) has some space-time-dependence, arising in part from the spatial variation 
of the background fluid temperature and velocity across the bubble wall, as discussed in 
Sec. 3. The background fluid is in local thermal equilibrium and comprises all ‘light’ effective 
degrees of freedom. Note that quasiparticles with large field-independent masses will be 
irrelevant for our purposes, since their distribution functions feature significant Boltzmann 
suppression. Also, precisely which fields should be considered ‘heavy’, ‘light’, or irrelevant 
depends on the given model. For the singlet-driven scenarios we are concerned with here, 
the heavy fields will be the top quarks, gauge, Higgs, and singlet bosons. 

To find approximate solutions to the Boltzmann equations for the heavy species and 
background, we will utilize the ‘fluid ansatz’ [77], in which case the perturbations are 
assumed to take the form 

~ ^( ST i + 5 T bg) ~ Pz(Svj + Vb g ). (4.6) 

Here fij. 6Tj, 5vj are the chemical potential, temperature perturbation, and velocity per¬ 
turbation of the species j, respectively, in the plasma frame. We have assumed that the 
fields with small couplings to the scalar condensates <fh,s are hi thermal equilibrium at a 
common space-time-dependent temperature T + +hTb g (x) and velocity fb g (x) with vanish¬ 
ing chemical potential, as in Ref. [77]. The assumption that /Ubg ~ 0 is valid whenever the 
total background particle destruction rate is larger than that for the heavy particles, as 
will be the case here (all pure gluon rates are enhanced by the large color factors and Bose 
statistics). The space-time-dependence in A7b g , Ub g arises from the change in masses of 
the corresponding particles moving from the <fi f 0 phase inside the bubble to the (j>i = 0 
vacuum outside. 

Throughout this study, we will work to linear order in the perturbations, which are 
assumed to be small (/ JLj/T , 5Tj/T , 6T^ g /T, 5vj, v^g -C 1). This should be the case 
for moderately strong phase transitions, and we verify the validity of this assumption a 
posteriori. It should be noted that this treatment can be extended to accommodate large 
fluid velocities in front of the wall [54], although this will not be necessary for any of the 
transitions we consider. As a result, we set all Lorentz 7 factors to 1 throughout our 
calculation. 

With the above definitions, the population factor V is given to linear order in the 
perturbations by 


'P — /o,l/o,2(1 ± /o, 3 )(l ± /o, 4 ) (< 5 l + $2 - £3 - 64) ( 4 . 7 ) 

where the ‘0’ subscript indicates the corresponding equilibrium distribution function. Note 
that the background temperature and velocity perturbations do not enter the collision 
integrals to linear order. 

To determine p,,;, <JT), and 5vi we follow Refs. [76-78] and take three moments of each 
equation, multiplying by f d 3 p/( 27r) 3 , f E p /Td 3 p/(2ir ) 3 , f p z /Tdfp/flirf and solve the 
resulting expressions for the perturbations. For a given heavy species, the relevant three 
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equations are given in the plasma frame by 


4* + 4 (Sr ' + ^ ) + f + « + j ( 2 jt) 3 T 2 


d 3 P C [f}. = c i dm i 


4^ + 4t m + ST ^ + CJ rm 


{SVi + U bg ) + J ~ 


2 T 9 t 
c| 9m? 
2T~df 


(4.8) 


f + + 5Tbg) + + Vbg) + / (£p|3 C[/]i 0 

where an ingoing particle of the relevant species has momentum p and where 


c l = 


g‘~ 2 f , d 3 p 

J'n+l 1 •'O ,i) ' 


(4.9) 


Further details can be found in Ref. [77]. The resulting collision terms for each heavy field 
i can be written as 

d 3 p 


2 (i w r‘„ ld + ST&.J 


d 3 p 


J^TsEiCWi ^E + STi^nj) 


d 3 p 




(4.10) 


The background excitations also satisfy a set of Boltzmann equations, 

Ed 3 p 


V- / 9 c 4 T 9 

C4 ( ai 6Tb z + 44dz Vhs ) + 




(2tt) 3 T 3 

Pzd 3 p 

(2vr) 3 T 3 


C't/lbg = 0 
C't/lbg = 0 


(4.11) 


which arise from Eq. 4.8 with ~ 0. The sum above is over all background species, with 
C 4 = J2 c 4 the heat capacity of the plasma. As for the heavy quasiparticles, the collision 
terms can be written as 


d 3 p 


( 2ir) 3 T 3 
d 3 p 


EiC[f] b g = - E + ^ r T 2 J 


(2vr) 3 T 4 


Pa.iC't/Jbg = - E (^Ai, 


(4.12) 


Although <5Tb g and Ub g do not enter the collision integrals, the perturbations corresponding 
to the heavy excitations do. The convention for evaluating the matrix elements is the 
same as for the heavy particles, with all background excitations treated as one species. 
Thus, every heavy particle process involving the background excitations will contribute to 
Eqs. 4.12. We will calculate all of the contributions relevant for singlet-driven transitions 
in the next subsection. 
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4.2 Relevant Excitations and Interaction Rates 


In the SM and its singlet extensions, the relevant heavy species to consider in Eq. 3.9 above 
are typically the top quarks, SU(2)l gauge, Higgs, and singlet bosons, with the Higgs and 
singlet excitations being the dominant source of friction on the singlet held condensate. 
We will consider two different sets of contributions to the total friction. For the gauge- 
invariant calculation, we include only the top quark, Higgs, and singlet contributions. When 
incorporating the gauge boson cubic term, we will also account for the friction arising from 
the SU(2)l gauge bosons. We will not include the Goldstone friction contribution, since 
we drop the corresponding finite-T cubic term from the effective potential. This should be 
a reasonable approximation as the Higgs excitations will only make up about 20% of the 
total friction on the wall 15 . 

The friction from each species enters the bubble wall EOM (Eq. 3.9) through the 
derivative of the corresponding mass squared. For the top quarks, the effective mass squared 
is 

+ n*(T) (4T3) 

with the corresponding thermal self-energy correction 

Ht(T) ~ l -glT\ (4.14) 


neglecting the subdominant thermal SU(2)l and U{1)y contributions. 

The Higgs and singlet require slightly more care. Throughout the remainder of this 
study we will neglect all mixing effects between the SM-like Higgs and singlet excitations. 
As discussed in Sec. 6.1, we will choose the parameters of the T = 0 Lagrangian such that 
the mixing vanishes in the broken electroweak phase. At finite temperature and across the 
bubble wall this will no longer be the case. However, in the high temperature limit, the 
effective neutral scalar mass matrix is diagonal, since off-diagonal thermal corrections are 
proportional to dimensionful parameters and vanish as T —> oo. For temperatures around 
the electroweak phase transition, the thermal masses still dominate the mixing matrix, and 
so this should be a decent approximation. The relevant field-dependent masses, including 
the leading thermal corrections, are then 


m 2 h (cl) h , (j> s ) ^ + 3A$ + -ai<j> 8 + ^a 2 4> 2 s + H h (T) 

m 2 s {(t>h, < h ) ^ b 2 + 26 3 </> s + 3 + n s (T) 

with the thermal masses 


(4.15) 


n h (T) ~ 
n s(T) = 


(^2 2 + ^i 2 + ^ 2 + ^A + ^a 2 ) 

G° 2 + i 64 ) t ' 2 


rp 2 


(4.16) 


15 We do include the Goldstones, c/> 0 , in the various interaction rates. They only appear as external 
legs of the diagrams we consider, and the corresponding matrix elements are gauge-independent. They are 
treated as a background species. 
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where we have left out the light fermion Yukawa contributions. Since we are neglecting 
the finite-temperature tadpole contribution to the effective potential we also drop the a\ 
terms in the masses above. This is required for consistency, since these terms are precisely 
those that give rise to the finite-temperature tadpole. Finally, as in Ref. [77], we treat the 
transverse SU{2)l gauge bosons as a single species W with field-dependent mass squared 

mwiM = \gWh • ( 4 - 17 ) 

Transverse excitations do not acquire a thermal mass at leading order in the couplings. 
Longitudinal modes obtain an effective thermal (Debye) mass at leading order, correspond¬ 
ing to the inverse screening length of electric potentials in the plasma [120]. This is given 
by 

m 2 D ,w(T) - (4.18) 

in the Standard Model. Since the gauge boson friction is dominated by very infrared 
excitations, only the transverse contributions will be relevant. 

Our strategies for dealing with each of these types of excitations will differ. As we 
will see below, the top quark and Higgs interaction rates are typically sizable, and so the 
collision term plays an important role in the corresponding Boltzmann equations. This is 
not expected to be the case for singlet quasiparticles at high temperature. Contrary to 
the tops and Higgs, we will assume that the singlet interactions are slow. In this case, the 
collision term can be neglected. The corresponding Boltzmann equation decouples from the 
rest of the system and can be solved exactly. We discuss this further in Secs. 4.2.2 and 5.1. 
Finally, the gauge boson contributions are dominated by infrared dynamics and require a 
classical treatment, which has been worked out in Ref. [77] and discussed in Sec. 4.3 below. 

Let us first consider the interactions involving the top quark, Higgs, and background 
excitations. 

4.2.1 Top, Higgs, and Background Excitations 

Solving the Boltzmann equations for the perturbations ST t <57b g , 5vt t h, and Ub g 

requires computing the collision integrals corresponding to all the four-body interactions 
involving t, h, and the background fields. This task is rather daunting due to the sheer 
number of allowed processes. However, the dominant interactions will be of 0(a 2 ) for the 
top quarks, and 0(a s at ), 0(a 2 ) for the Higgs bosons, where a s = g 2 /An, at = y 2 /47 T. We 
will therefore focus on these interactions, neglecting, for example, contributions involving a 
factor of a w , which are numerically small compared to the Yukawa-type contributions for 
the Higgs bosons 16 . 

To estimate the relevant interaction rates, we will work at leading order in all couplings 
in the high-T, weak coupling limit, neglecting all terms of 0{m 2 /T 2 ) (here m should be 
understood as either a zero-temperature or thermal mass). This is the approximation 
used in all previous microphysical studies of the wall velocity [76-78], as well as in the 
context of plasma properties in arbitrary high-temperature gauge theories [121, 122], This 

1G We have verified this is the case despite the enhancement provided by Bose-Einstein statistics. 
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approximation can begin to break down inside the bubble wall for the top quarks and scalars 
and could in principle be improved upon in the future. Nevertheless, it should reproduce 
at least the correct parametric dependence of the full leading order result. In this limit 
we can neglect the effect of the space-time-varying masses on the interaction rates. Hard 
excitations dominate the phase space for the relevant top quark and Higgs collision integrals, 
and can be characterized by massless dispersion relations 

£ ,hard = yjp 1 2 + m 2 + m 2 h (T) ~ p + O ( gT ) (4.19) 

to leading order in the couplings and in the high-temperature limit. 

Although the external quasiparticles can be treated as massless in this approximation, 
infrared excitations appearing as mediators in t- and u-channel diagrams naively result 
in logarithmic IR divergences in the Higgs and top quark scattering amplitudes. These 
divergences are cut off by the interactions of the mediator with the plasma. For long- 
wavelength excitations, the corrections comprise so-called ‘hard thermal loops’ (HTLs), and 
result in a breakdown of the perturbative expansion. The corrections can be resummed into 
a thermal self-energy correction to the propagator, valid in the low momentum limit 17 . The 
self-energy is typically of order gT and so these processes can produce sizable logarithmic 
enhancements of the corresponding matrix elements, scaling as g{T) 4 log l/g(T) at high 
temperatures. This provides a useful way of categorizing the most important diagrams 
contributing to A4* in C[f] in the high-T limit. 

A full leading-order determination of the effective scattering rates in the plasma is 
possible [122, 123], though computationally more involved and beyond the scope of this 
work. It would be interesting to revisit in the future. We will instead work in a ‘lead¬ 
ing logarithm’ expansion, keeping only contributions of order rs_/ 5 (T) 4 loglMT), which 
are typically the largest. In this approximation, only 4-body rates with t- and tt-channel 
diagrams contribute. For further details on this approximation, see Refs. [77, 121], 

Another subtlety arises in computing scattering rates in the high-T limit involving 
soft t- or it-channel exchange. The thermal self-energies involved in the propagators are 
generally momentum-dependent [122], Previous studies of the wall velocity neglected these 
contributions, simply replacing them with the corresponding Debye (thermal) masses for 
the corresponding gauge bosons (fermions). However, including the momentum-dependent 
self-energies, which enter at leading order in the couplings, can have a significant effect 
and has been shown to often provide better agreement between the leading log and full 
leading order results for plasma transport coefficients in high-temperature gauge theories 
[122], Consequently, we will use the full momentum-dependent HTL-resummed propagators 
[122, 124, 125] in computing the various collision integrals for the top quark and Higgs 
excitations. 

The relevant processes and their associated vacuum matrix elements are listed in Ta¬ 
ble 1. All terms of 0(m 2 /T 2 ) have been dropped. The leading log matrix elements are 
summed over the helicities and colors of the external particle (but not particle-antiparticle). 

1 'Note that the same situation arises when computing loop corrections to the finite-temperature effective 

potential. 
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Process 

\m\L 

Internal Propagator 

O(gi): 
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tg 77 tg : 
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160 gl^ 2 

9 

o(vhlY 

tt 77 hg, (jpg: 
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t 
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(¥ + B 

t, b 

tg 77 th , 

-SVtffii 

t 

tg 77 bcj) + : 


b 

t(j)~ 77 bg: 

~*Vt9 If 

t 

0(ytY 

tt 77 hh, (j)°(j) 0 : 

bt (7 + i) 

t 

tt 77 (f) + (j)~: 

zyfi 

b 

tt 77 h(f>°: 

bt (7 + h) 

t 

tb 77 h(/) + , 4>°4> + - 

bti 

t 
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_3 „ .4 s 

2 yt t 

t 

tcj>~ 77 hb , 4>°b: 

bti 

t 

tct) + 77 4 >+f: 

3% 4 7 

b 


Table 1. Relevant 4-body processes and their corresponding matrix elements in the leading log 
approximation. The matrix elements are summed over the helicities and colors of all four external 
states (as well as flavors and quark anti-quark for tq —y tq ). The excitation appearing on the 
internal propagators in each case is listed in the right-hand column. Note that other t- and u-channel 
processes exist, but do not contribute logarithmically to the collision integrals or are suppressed by 
powers of couplings small compared to 53 , yt- 


These contributions will also be divided by the number of degrees of freedom of the species 
under consideration when entering the various Tji, ■. 

The vacuum matrix elements must be modified to include the medium-dependent effects 
discussed above. At leading order, this amounts to inserting the momentum-dependent HTL 
self-energies on the internal lines. To translate the vacuum matrix elements above to their 
finite-temperature analogs, we can use the results of Refs. [119, 122], For diagrams with an 
exchanged fermion in the leading log approximation, this amounts to the replacement 

u _ _s 4:Re{p-gk-g* +sq-g*) 
t~ t | q-q\ 2 
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with = p^ — p ,fl — £^(p — p r ) and E^g) the fermionic HTL self-energy function [122, 
124, 125] 


S°(?) = 
S(9) = 


m?AT) 

2|q| 



q 


|q| + q° 

|q| - Q° 

1 + in — 



2|q| 


log 


q| + q° \ 
q| -q°) 


(4.21) 


with rrif(T) the leading order fermion thermal mass, given approximately by rrif(T) ~ 
yJ\/Qg$T for quarks (see Eq. 4.16). 

For the gluon exchange diagrams, we must replace 


s 2 + u 2 1 


t 2 


-(i+\D^(p-p')( P +p'nk+kr\ 2 ) 


(4.22) 


with the retarded thermal equilibrium gluon propagator D^^q) given by 


D oo(q) - 2 

|q| +n 0 o(g,r) 

D (a) = 5ij ~ 

g 2 + n r (v,T) 

Ao(<?) = D i0 (q) = 0. 


(4.23) 


The relevant HTL gauge boson self energy is [122, 124, 125] 


(4.24) 


with mn{T ) the Debye mass of the gauge boson (= \[2g^T for the gluon). 

With the above replacements, we can now perform the collision integrals numerically. 
We do so using the phase space parametrization discussed and detailed in Refs. [119, 122, 
126] and the Vegas Monte Carlo routine included in the Cuba package [127]. The integrals 
converge reasonably quickly for most cases on a standard desktop computer to reasonable 
numerical precision. 

The results of this numerical evaluation are given in Equations 4.25-4.29 below. These 
values can then be plugged into Eqs. 4.8 and 4.11 for the perturbations. For the Higgs 
bosons, the rates (in the leading log approximation) are 


n 0 o(<?) = rn^(T) ( 1 - A log -% + 


2 q 


nr q 

2]qj 


I Wq)=m&T) 


+ 


A 2 


2|q| 4 |q 


i q +<f 

log rn-n - m 


T ^ (El X 10 ~ 3 g 2 y 2 + 6.0 x 10 ~ 4 yf)T 
^Ti,h - r ^2 ,h - ( 2-5 X 10 ~ 3 g 2 y 2 + 1.4 x 10 ~ 3 y 4 )T 
r^, fc - (8.6 x 10~ 3 g 2 y 2 + 4.8 x W~ 3 yt)T 
r- (3.5 x 10~ 3 g 2 3 y 2 + 1.8 X 10~ 3 yf)T, 
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while the corresponding contributions to the top quark distributions are 


- rj 1|t ~ (1.0 x 10- 3 <7 ly 2 t + 5.8 x 10 ~ 4 yf)T 
~ r Ti,t - r^ 2 ,t - (2.5 X 10+ 1.5 x 10 ~ 3 yf)T 

- V^ 2 t ~ (8.5 x HT 3 ^ 2 + 4.8 X W~ 3 y?)T 

- rj )t ~ (2.8 x 10 ~ 3 glit + 1.4 x 10~ 3 yf)T. 

For the top quarks, 

r^ 1>t ~ (5.0 x 10- 4 ff | + 5.8 x l(T 4 5 2 y 2 + 1.5 x Kr 4 y 4 )T 
rKt - ,t - (1-2 X 1(T 3 5 4 + 1.4 x 10 ~ 3 g$tf + 3.6 x 10"V)T 

^T 2 ,t - (1-1 x 10- 2 5 4 + 4.6 x 10- 3 5 |y 2 + 1.1 x l(T 3 y 4 )T 
Ft lt - (2.0 x l(T 2 5 f + 1.7 x 10- 3 ff 2 y 2 + 4.3 x l(r\ 4 )T, 

while their contributions to the Higgs distributions are 

- T U,h - ( 9 - 3 x w ~ 5 dVt + 5-3 x 10~ 5 y 4 )T 

- r Ti,h - r U,h - (2-2 x 10- 4 g 2 y 2 + 1.3 X 10-\ 4 )T 

- T l T2th ~ (7.2 x 10 ~ A gly 2 t + 4.0 X 10-\ 4 )T 

- Y^ h ~ (2.4 x 10- 4 5 2 y 2 + 1.2 x l(T 4 y 4 )T. 

Finally, the background contributions are 

,t - (1.4 X 10 _ 2 g| + 1.3 x 10~ 2 g|y 2 + 2.6 x 10~ 3 y?)T 
f T 2 ,t ^ (1.4 X 10 —1 <73 + 4.6 x 10- 2 ^ + 8.7 x l(T 3 y 4 )T 
r v ,t ^ (2.4 x 10 _ 1 g| + 1.7 x 10- 2 g 2 y 2 + 3.4 x l(T 3 y*)T 

r M 2 ,h - o 

r T 2 ,h - (1.0 x 10 _3 g|y 2 + 9.8 X 10 _5 y 4 )T 
f v>h - (1-6 x 10- 3 .g|y 2 + 4.6 x l(rV)T. 

Note that the background rates for the top quarks tend to appear larger than their counter¬ 
parts in Eq. 4.27 above. This is simply because in the background rates we have summed 
over all contributions, while the rates in Eqs. 4.25-4.28 are the average values per degree 
of freedom (e.g. divided by Nt = 6 in the top quark case). The latter rates will be mul¬ 
tiplied by the appropriate iVj factors when they enter the bubble wall equation of motion. 
Also, note that the contributions in Eqs. 4.26 and 4.28 are negative because they arise from 
diagrams with the relevant species on the outgoing legs of the Feynman diagrams. 

In previous work, the above integrals were performed analytically using several approxi¬ 
mations and without incorporating the (momentum-dependent) self-energies. The different 
computational methods used here change the rates by 0(1) factors relative to the results in 
Ref. [77]. There were also some algebraic errors in the results of Ref. [77], as pointed out in 
Ref. [121], that contribute to the discrepancy. Although our treatment is still formally at 
the same order as that of Ref. [77], the HTL-improved calculation in many cases is expected 


(4.26) 


(4.27) 


(4.28) 


(4.29) 
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to more closely reproduce the leading order result (see e.g. Fig. 1 of Ref. [122]). Neverthe¬ 
less, modulo the algebraic mistakes in Ref. [77], our interaction rates are no more accurate 
than those of Moore and Prokopec in approximating the full leading order results; they 
should simply be thought of as arising from a different set of approximations. We comment 
further on the differences between the rates found in Ref. [77] and those reported above in 
Appendix A. The reader should bear in mind that the predicted wall velocity will tend to 
be higher if the rates computed in Ref. [77] are used instead of ours. This is because the 
former are larger and thus result in faster equilibration for the various perturbations. Note 
also that the collision integrals computed and listed above depend only on the Standard 
Model degrees of freedom, and as such are quite general. They can be used in various 
applications beyond those considered in this work. 

Before moving on, some comments regarding the higher order contributions neglected 
in the fluid approximation are in order. The assumed form for the perturbations is that of 
a perfect fluid and can be thought of as a truncated expansion in powers of momentum. 
That is, the fluid ansatz assumes that the effects of higher angular moments p e Y( m { p) hi 
the distribution functions are negligible [77] (the Y( m are spherical harmonics). For the 
top quarks this is a good approximation, since we find that the velocity perturbations 
typically satisfy dvT/dp < 0.1, while the contributions from higher moments should scale 
roughly like (6vT / 6fif [77]. On the other hand, the Higgs bosons have smaller interaction 
rates than the tops, and so the fluid approximation begins to break down for strong phase 
transitions. For this reason we will restrict ourselves to moderately strong phase transitions 
with <ph{T n )/T n < 1.1 in our consideration of the xSM in section 6. Already in this regime 
some points will be found to possess no subsonic solutions. Further details on the limitations 
of the fluid approximation can be found in Appendix B of Ref. [77]. 

4.2.2 Singlet Contributions 

Excitations of the singlet held will also contribute to the friction on the bubble wall. The 
corresponding collision integral for singlet quasiparticles is dominated by scattering pro¬ 
cesses involving four external scalars. At high temperatures, the resulting effective inter¬ 
action rates are typically suppressed relative to those for the processes involving external 
fermions. To see this, note that processes with f-channel diagrams involving two external 
scalars and two external fermions schematically contribute 



to the first moment of the Boltzmann equation in the small-g limit (here q = |p — p'l). The 
logarithmic divergence is cut off by the thermal self-energy of the exchanged quasiparticle. 
The upper limit q ~ T corresponds to the breakdown of the small q approximation. For pro¬ 
cesses involving four external scalars with a scalar exchanged in the f-channel, the integrals 
of the Bose-Einstein distribution functions are also IR sensitive. Cutting off the distribution 
functions with a parameter e with mass dimension 1 such that fo(p/T) —> fo{p/T + e/T), 
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we find the corresponding contribution to be 

a 4 T f T , 1 . 

r M .i ~ —5- x / (4.31) 

e 2 im(T) r 

in the small q and e/T regime (a is a cubic coupling with mass dimension 1). The diver¬ 
gence is now nominally quadratic (again cut off by the self-energy of the exchanged scalar), 
and the integral of over the distribution functions is cut off by the thermal masses of the 
external scalars in the infrared. This suggests that a rough leading order estimate of the 
scalar quasiparticle scattering rates should include thermal masses in the Bose-Einstein 
distribution functions. Computing the dominant contributions involving the cubic and 
quartic couplings and performing the resulting integrals, we find interaction rates that are 
significantly smaller than those for the tops and Higgs across the range of couplings and 
temperatures we consider, despite the nominally more severe divergence structure. This is 
because for large temperatures, e ~ gT and so the schematic rate in Eq. 4.31 is suppressed 
by 1/T 3 . This is expected, since at high temperatures all dimensionful parameters of the 
zero temperature theory should be irrelevant [123]. Meanwhile the quartic interactions 
do not contribute a small q divergence at leading order. We therefore expect the singlet 
qusiparticle collision term to be small. 

Assuming this is the case, the fluid approximation is likely to provide a rather poor 
estimate of the corresponding friction. Instead, we will make a ‘free particle’ approximation 
[72, 73] for these excitations, dropping the collision term for the singlet. In this case, the 
Boltzmann equation can be solved exactly, without taking moments. The solution is given 
by Eq. 5.3 of Ref. [77], to lowest order in v w , and is reproduced below in Eq. 5.1. We 
will include the corresponding contribution to the equations of motion when computing the 
friction. Note that the presence of non-negligible interactions would decrease the friction 
and increase the predicted wall velocity. 

This treatment assumes that the dominant singlet excitations are well described by our 
perturbative kinetic theory, with a local collision term. This should be true for the hard 
excitations. For much softer excitations, a classical treatment with the short wavelength 
fluctuations integrated out is likely more appropriate [79, 128]. Unlike classical Yang-Mills 
fields (discussed in the next subsection), the classical scalar field is not overdamped [129], 
and so infrared excitations are likely to equilibrate quickly. We thus neglect the effect of 
infrared singlet modes on the friction. This approximation is likely rather crude and should 
be revisited in the future. Including the IR contributions would increase the friction and 
potentially slow the wall down. The reader should bear this in mind as we proceed. 

4.3 Gauge Boson Contributions 

Finally, for our calculations incorporating the finite-temperature gauge boson cubic term 
in Eq. 2.7, we will include the friction from the SU{2)l gauge bosons. In contrast with the 
top quarks, Higgs, and singlet excitations, the friction in this case is dominated by infrared 
degrees of freedom, which can be treated approximately as over-damped classical fields [129] 
as opposed to the perturbative approach utilized for the other species. The distribution 
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functions in the classical limit can be shown to satisfy [129] 

™hMT)d M r,T) = _ [p , + mU<Ph)]MpX)+M ( 4 . 32 ) 

where J\f is a noise term. This equation serves as the analog of Eq. 4.2. Further discussion 
of its derivation and applicability can be found in Ref. [129]. 

Note that hard gauge boson excitations also exert a drag force on the bubble wall, as 
computed in Ref. [77]. However, we have verified that these contributions are substantially 
suppressed relative to that from the IR gauge bosons, as found in Ref. [129]. We do not 
include them. 


5 Solving for the Wall Velocity 

With the collision terms evaluated, we can now solve the Boltzmann equations to determine 
the perturbations for a given held profile and wall velocity. The goal is then to find the 
value of v w and the profile (and hence the perturbations) such that the equations of motion 
are satisfied. We will describe how this can be done below. First, let us consider solutions 
to the Boltzmann equations given a particular profile and value of v w . 


5.1 Exact Solution for the Singlet Excitations 


As mentioned above, the singlet equation can be decoupled from the rest of the system 
and solved exactly in the free-particle limit. The result been discussed in detail previously 
[72, 73, 77], and so we simply quote it here. The integral appearing in the equations of 
motion 3.9, to lowest order in v w , is 


f d 3 p f d 3 p e E p/ T Q(Pz) 

J (27r) 3 2E 5fs ^ X) ~ Vw j (2it) 3 2E ( e E p /T ± X ) 2 T 


(5.1) 


where the upper (lower) sign is for fermions (bosons). The function Q is defined as 


n , •, _ j Vp 2 z + m s ((t> h ,(j) s ,T) 2 -p z , p z > -^mP s (T) 2 - m s ((j) h , (f> s , T) 2 

{Pz> \ ~VpI + <l>s, T ) 2 - m2(T) 2 - p z , p z < - v / mO(T)2-m s (^,^,T)2 

(5.2) 

with mP s {T ) the singlet mass (including thermal contributions) in the broken phase (i.e. 
at z —> oo). This integral is multiplied by dm 2 s {(j)h, </>s,T) /dcj)h hi the Higgs EOM, and by 
dm 2 s {(f) h ,(j) s ,T)/d(j) s in the singlet equation. 


5.2 Exact Solution for the IR Gauge Contributions 

We can also solve for the classical gauge boson contribution in Eq. 4.32. The result is [129] 
drriyyi&h) f d 3 p 3T 2 (/)' h (x) 

J Wm ,w(p ' x) = 0(1 “ x<) (5 - 3) 


where the quantity x* solves mwifihix*)] = 1/h/i, with Lh the SM-like Higgs wall width. For 
smaller x, the WKB description used to derive Eq. 4.32 breaks down. For more discussion 
on this point, see Ref. [129]. Note that this value cuts off the IR divergence of Eq. 5.3. 
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5.3 Solving the top-Higgs System 


It remains to solve the equations for the top quark, Higgs, and background excitations. 
Here we follow the methods found in Refs. [77, 78] with some slight modifications. 

Since we are interested in static solutions to the equations of motion in the wall frame, 
all quantities depend only on x and so the derivatives in Eqs. 4.8 can be re-written as 
dt —> v w d/dx, d z —> d/dx. The Boltzmann equations in the static limit are therefore a set 
of linear ordinary inhomogeneous differential equations. 

To solve them, the equations for the background temperature and velocity, Eqs. 4.11, 
can be used to eliminate T^ g and Ub g from the top quark and Higgs equations. Defining a 
vector of perturbations 

5 t = (5/j.t, 5T t , 5v tl 5Hhi 5T h , 5v h ) , (5.4) 

Eqs. 4.8 can then be written as 


Ak~rh + Ti k 5k = Fi 
dx 


(5.5) 


with the definitions 
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(5.6) 


(5.7) 


and the source vector 


F(x) 1 = 


X _ 

2T 


dml((t> h ) 

dx 


dm ? 




dx 


0, cf(x) 


dm 2 h ((j) h , 4> s ) 
dx 


, dm\(cj)h, 4>s 
dx 


(5.8) 


The field-dependent masses are given by Eq. 4.15. 

The system of equations can be solved by simple Green’s function techniques. Following 
Ref. [77] we define the matrix x such that 


{A 1 r). j Xjk = Xikh (5-9) 

where X k are the eigenvalues of H _1 T. It is then straightforward to define the vector Green’s 
function 

Gi(x, y ) = sgn(Ai)e _Ai(a:_?/) 0 [sgn(A;)(x - y)\ (5.10) 

in terms of which the perturbation 5i is given by 

/ OO 

[x-'A-'Fiy)]. Gj(x, y)dy. (5.11) 

-OO 


27 - 









These solutions can then be inserted into the equation for the variation in the background 
temperature, 


5T hg (x) = / V TY v j8vj - v w (+ T T2) i5TA ( 5 . 12 ) 

c 4 (3 -vi) J-oc^ 1 V yj 


defined with respect to T+, its value far ahead of the bubble in the shock front. 


5.4 Approximate Solutions to the Equations of Motion 


With the perturbations determined, we can now try to identify solutions to the wall equa¬ 
tions of motion. In terms of the perturbations Sj, Eq. 3.9 reads, for the gauge-invariant 
case, 


,, 2 \ jji , dVfaT) , ^dm^i)T 

■(! “ v w)<Pi + —TZ—+ 2^ 


d(t>i 


T 


dcj)i 

dm 2 (fa) 
d(pi 


c{5pj + c^(6Tj + ST bg ) 

d 3 p 


(5.13) 


(2vr) 3 2 E 


$fs{x,p) = 0 


where the last term is given in Eq. 5.1. If the gauge boson contributions are included, 
the RHS of Eq. 5.3 should be added to the LHS of the above expression. The bound¬ 
ary conditions are 4>h,s(x —> Too) = 0ft, s; ±(T+) and 4>' hs {x —> ± 00 ) = 0. This system of 
equations will typically admit a solution for certain values of v w and profile (j>h,s(x)- Our 
strategy will be to vary the profile and scan over values of v w consistent with a deflagration 
bubble, looking for parameters such that the equations of motion (and Boltzmann equa¬ 
tions) are simultaneously satisfied. All parameter space points we consider have at most 
one deflagration solution. 

Eq. 5.13 represents a set of integro-differential equations for $ = ((fih, 4> S ) T (in the fol¬ 
lowing discussion it will be useful to use explicit vector notation). To find its approximate 
solutions, we will follow the strategy of Refs. [77, 78] and use an ansatz for the field profiles 
which depend on only a few parameters. Of course in using an ansatz it is unlikely that the 
full equations of motion will be satisfied exactly. However, we can reasonably approximate 
a solution by scanning over the ansatz parameters and imposing physical constraints. For a 
given choice of parameters, the Boltzmann equations can be solved exactly and the results 
inserted into the EOM. A set of parameter values such that all constraints are simulta¬ 
neously satisfied corresponds to an approximate solution to the original equations. This 
strategy has been employed in previous calculations of the wall velocity [54, 76, 77] and we 
expect the results obtained in this way to be a decent approximation to the full numerical 
solution. 

Before analyzing Eq. 5.13 further, some useful insight can be gained from solving the 
corresponding field equations with the assumption of constant friction of the form 

„ 4 2 a> 4 * 

+ + ^ = 0 < 5 ' 14 ) 

subject to the boundary conditions <3?j(x —> T°°) = d>j.±, 4)'(x —> Too) = 0, where 
= ( d/d(j)h,d/d(j) s ) rT and T is the same for both field directions. Clearly this is not 
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a realistic case, but we will improve on it below. These equations are much simpler than 
the full integro-differential equations of Eq. 5.13. As detailed in Ref. [48], the solution to 
the equations of motion can be found numerically via path deformations, and corresponds 
to a limit in which all the friction on the wall is parallel to the trajectory in the (fih, fi s ) 
field space within the wall. This can be seen by noting that d<&/dx is a “velocity” vector in 
field space, so J-d^/dx always acts parallel to the trajectory. In all cases we consider, the 
solution to these equations of motion is well-fit by a tanh ansatz with parameters 

fii(x) = -y + tanh X (5.15) 


in the fluid frame. Here Li = Lf ,, L s are the wall widths and Si are offsets allowing for a 
good fit to the numerical solution. We can take <5/,, = 0 without loss of generality. Note 
that if we had allowed (f> s 0 in the electroweak-symmetric phase, we could have instead 
used (j) s {x) = + A<f s /2 ^1 + tanh for the ansatz, with 0° the singlet VEV in the 

symmetric phase and Acf> s the change in VEV across the wall. The remaining analysis 
would proceed in the same way. 

How does the situation change when including a realistic friction term? The friction 
is no longer proportional to d^/dx and so will have some component perpendicular to 
the field space trajectory, acting as an effective “normal force" along the path. However, 
there is a fortunate simplification we can make if the friction perpendicular to the held 
space trajectory found by solving Eq. 5.14 is negligible. In this case, the field will not be 
significantly deformed from its field space path found using the constant friction equations of 
motion , although the held prohle in physical space will change. In other words, if we write 
the solution to Eq. 5.14 as < J?(s) where s = s(x) is some parameter such that \d&/ds\ = 1, 
the effect of a change in the friction parallel to 3?(s) will only be to alter s(x). Meanwhile, a 
change in the friction normal to the prohle would result in a change of ^(s) itself. Applying 
this reasoning to the tanh ansatz (which we find to be a good ht to the constant friction 
solution), the effect of altering only the friction parallel to the trajectory and neglecting 
that normal to the path will simply be an overall simultaneous re-scaling of all the widths 
and offsets: 

Lh,s t aLfo^si S s y aS s . (5.16) 

This is the only change in the tanh prohle that will not deform the path in held space. 
Then, starting from the constant friction solution, the problem can be reduced to finding 
the values of v w and a such that the pressure and pressure gradient in the wall vanish: 


d& 

dx (Eq. 5.13) • — = 0, 
dx 


J dx (Eq. 5.13) 


d 2 $ 

dx 2 


= 0 . 


(5.17) 


The above constraints will only be satished for values of v w and a such that the wall is 
not accelerating or expanding/contracting, as required for the steady-state solution we are 
seeking. This is a simple generalization of the strategy used in the SM case in Refs. [54, 
76, 77], where v w and L w are varied. 
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In what follows, we will assume that the friction force normal to the field space path 
determined from Eq. 5.14 is negligible, such that the discussion of the above paragraph 
applies. The validity of this assumption can be checked a posteriori (which we do), but 
there is an intuitive reason why it should often be reasonable. At very high temperatures the 
potential is stabilized at the origin (in our approximation) by the effective thermal masses 
~ gT of all the scalars. Around T ~ 0, the potential is necessarily stabilized at a minimum 
away from the origin. The temperature of the phase transition is such that two minima 
exist simultaneously and are (nearly) degenerate. Provided that the tree-level contribution 
to the barrier is not too large, this can only occur if there is a significant cancellation 
between the zero-temperature and finite-temperature corrections in some direction of the 
held space. This approximate cancellation will be largest along the held space trajectory 
<&(s) found by solving Eq. 5.14 such that, schematically, 


v^(*,r = o) 


d&(s) 

ds 


E 


dm(<&Y 

d&(s) 


d 3 p 

( 2ir) 3 2E 


fo (P,T) 


d$(s) 

ds 


(5.18) 


The resulting ridge in the hnite-temperature effective potential is precisely that along which 
the cubic term becomes relevant. We can insert the solution to Eq. 5.14 into Eq. 5.13 
and see how we expect the solution to change when going to the full EOM. With the 
approximate cancellation of Eq. 5.18 in effect, the full equation of motion parallel to <&(s) 
is then schematically 


. o , d? $ v—' 

~^ 1 ~ v ^lx 2+ ^ 


dm(i f>) 

d$ 


d 3 p 


(: 2ir) 3 2E 


5f(p,T) + 0(5f 2 


d${s) 

ds 


(5.19) 


along <&(s). The (approximate) cancellation has made the contribution from the friction 
the leading effect. The change in the friction term from Eq. 5.14 to 5.13 will alter s(x) but 
leave the field space trajectory <&(s) unchanged. In contrast, the cancellation in Eq. 5.18 
is not expected to hold in the perpendicular direction along the path (oc d 2 &/ds 2 ), unless 
the minimum of the potential away from the origin lies in the bottom of a shallow valley. 
In the absence of such an approximate continuous symmetry, the full EOM perpendicular 
to 3>(s) is 


9 .d 2< ]> _ dm(<&) 2 

“( 1 ~ Vw ^~dx2 + = °) + 


d 3 p 


d$ 


(27t) 3 2£ 


f 0 (p,T) + O(6f) 


d 2 $(s 

ds 2 


(5.20) 

and so the effect of the friction in this direction is perturbatively small, resulting in a 
negligible correction to the perpendicular component of Eq. 5.14 and hence to <l>(s). 

The above discussion suggests the following strategy for Ending approximate solutions 
to the equations of motion for a given parameter space point: 


1. Compute the phase transition properties, namely the order parameter and T n . We 
do this using the CosmoTransitions package [130]. 

2. Solve for the constant friction profile from Eq. 5.14. This can be done using path 
deformations [48] or otherwise. Fit the solution to the tanh ansatz Eq. 5.15. 
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3. Solve the hydrodynamic relations to obtain T + for various values of v w . 

4. Vary the values of v w and a. For each pair, solve the Boltzmann equations as discussed 
above, using T = T + and Li = aL®, 5 S = ad|? with L? and <5j? obtained from the 
numerical solution of Eq. 5.14. 

5. Insert the solutions for the perturbations (and background temperature profile) into 
Eq. 5.13, then compute the constraints in Eq. 5.17. The values of v w and a satisfying 
Eq. 5.17 can be found by interpolating between the results of the scan. 


This method generalizes that of Refs. [54, 76, 77] to accommodate the additional singlet 
held direction. 

The results for v w and a obtained in this way will still produce a residual ‘normal force’ 
perpendicular to the trajectory in held space when inserted back into Eq. 5.13 due to the 
neglect of the friction in the direction oc dP^/ds 2 . Defining s(x) = |<I>(x)|, the tangent and 
normal unit vectors to the held space path <&(s) are 


t(s) = 


d&(s) 

ds 


d&(s) 


ds 


-l 


n(s) = 

ds 


(5.21) 


The ‘normal force’ along $(s) is given by 


d 2 <& ( ds\ 

N(i)= ^u) - 


v_w(4>,r) + ^r 


dm(<&)' 


d 3 p 


(! 2ir) 3 2E 


Sf(p,T) 


• n[s) 


(5.22) 


A full solution to the equations of motion would guarantee that N(x) = 0 for all x. This 
will not be true for our approximate solutions. 

To check that this residual normal force is indeed negligible, we can deform the prohle 
to eliminate it. The deformation can be performed along the lines suggested by Ref. [130] 
for computing the critical bubble prohle. It typically results in small changes to the original 
held space prohle, which in turn have very little effect on the perturbations and constraints 
in Eq. 5.17, since the curvature perpendicular to the path is typically significant. This 
suggests that the wall velocity and prohle found in the way outlined above should indeed 
provide a reasonable approximation to those obtained from the full solution of the equations 
of motion, at least in the cases we consider. There may be exceptions elsewhere in the 
parameter space. 

Note that the procedure outlined in Steps 1-5 above is quite general, and can be adapted 
beyond singlet models to other scenarios with multiple held directions, provided Eq. 5.18 
is approximately satished. 


6 Wall Velocities in the Real Singlet Extension 
6.1 Parameter Space and Phenomenology 

We now turn to the parameter space of the real singlet extension of the Standard Model as 
an application. Our goal is not a comprehensive analysis of this model. Instead, we focus 
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on a sample of the parameter space consistent with current observations and a strongly 
first-order phase transition. 

The authors of Ref. [49] recently performed a detailed analysis of the electroweak phase 
transition in this setup and so we use their study as a guide. Recall that we have identified 
the excitation h as the Standard Model-like Higgs with rrih — 125 GeV. The couplings of 
the discovered Higgs are very close to those expected in the Standard Model [131, 132], 
Therefore we will assume that the doublet H couples precisely as in the Standard Model 
and that there is no mixing between the singlet and Higgs at tree-level 18 . This corresponds 
to the choice 

a\ + 2 a 2 v s = 0 (6.1) 

and immediately fixes A in terms of the observed Higgs mass, m 2 = 2Xv 2 = (125 GeV) 2 . 
Note that Ref. [49] showed that deviations from this no-mixing limit are allowed by current 
LHC measurements and so this requirement can be relaxed. In our scans, we will vary both 
the cross-quartic coupling, a 2 , and the zero-temperature singlet VEV, v s . Then Eq. 6.1 can 
be used to determine a i, given our choice for v s and a 2 . 

We will also assume 63 = 0. The barrier required for a first-order electroweak phase 
transition will then arise primarily from the tree-level mixed cubic coupling a\. Again, this 
is not required by the phenomenology, and this choice can be altered without significantly 
affecting any of our arguments. Note that much of the NMSSM parameter space compatible 
with a strongly first-order electroweak phase transition lies close to the corresponding limit 
\kA k \ <C |AA a | [47, 48], 

Equating the minima of the tree-level potential with the VEVs v = 246 GeV and v s 
yields the conditions 

/u 2 = Xv 2 + -(01 + a 2 v s )v s 

n v 2 1 (6 ' 2) 

62 = -b 3 v s - b 4 v 2 - -j- - -a 2 v 2 

4v s 2 

which allows us to solve for ju 2 and b 2 . We also take the tree-level singlet mass squared (in 
the zero mixing limit), 

m 2 s = b 3 v s + 2b 4 v 2 - (6.3) 

as an input, and use the above expression to solve for 64 . 

The free parameters are thus a 2 , v s , and m 2 . There are some additional requirements on 
the theory that allow us to hone in on phenomenologically viable values for these quantities. 
First of all, stability of the T = 0 potential requires 64 > 0, which in turn limits the values 
of rn 2 we can consider via Eq. 6.3. Also, if m 2 is too small, the decay h —> ss would cause 
large deviations in the width of h which are not observed experimentally. On the other 
hand, if m s > 2m.h, di-Higgs production at the LHC can place significant constraints on 
the model [133, 134], For simplicity, we avoid this regime and choose mh/2 < m s < 2rrih- 

Further insight can be gained from considering the expected behavior of the electroweak 
phase transition strength and wall velocity as a function of these free parameters. In 

18 Departing from this choice should not significantly affect our predicted range of wall velocities, since 
we drop the finite-temperature tadpole and thus the contribution from 01. 
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particular, larger values of 02 are favorable from the standpoint of small wall velocities. 
This is because increasing 02 corresponds to lowering T c , as can be seen by noting that 
T c corresponds approximately to the temperature such that m|(0/j, <j> s ,T) ~ 0. Lowering 
T c results in lower values of <f> c satisfying 4> C /T C > 1. Smaller field values in turn lower 
the total pressure difference between the vacua. This pressure difference is what drives the 
expansion of the bubble and which must be compensated for by the friction. Furthermore, 
larger values for 02 correspond to larger m^((f)h, 4>s)i which in turn increases the friction 
from the singlet on the wall. 

With the above reasoning in mind, we choose two sets of parameters likely to be 
promising for electroweak baryogenesis and across which we can compute the bubble wall 
velocity. These are 


Set 1 : m s = 170 GeV, a 2 = 0.9 
Set 2 : m s = 245 GeV, a 2 = 1.7. 


(6.4) 


For both sets of points we vary v s , which corresponds to varying the strength of the elec¬ 
troweak phase transition. This is clear, since higher values of v s correspond to larger |ai| 
via Eq. 6.1 and hence a larger contribution to the barrier separating the electroweak min¬ 
imum from the origin at finite temperature. We vary v s up to values such that either the 
fluid approximation breaks down, or subsonic solutions no longer exist. This corresponds 
roughly to values of v s between 30-100 GeV for Sets 1 and 2. Note that with our choices of 
parameters, every point in Sets 1 and 2 will satisfy all current phenomenological constraints. 
Specifically, the electroweak vacuum is absolutely stable for all points considered while the 
absence of s — h mixing ensures that both h and the new singlet-like state are compatible 
with current observations and limits. 


6.2 Results 

Finally we arrive at our results for the xSM. The wall velocities computed in a parameter 
scan for the two sets of points (Set 1 and 2) described above are shown in Fig. 2. The critical 
bubble profile and nucleation temperature are computed using CosmoTransitions [130]. 

The solid lines depict the outcome of the gauge-invariant method. Stronger transitions 
correspond to faster moving bubble walls. The perturbative fluid approximation becomes 
worse and breaks down for stronger phase transitions, and so we cut off our scans above 
4>h{T n )/T n ~ 1.1. Wall velocities for stronger phase transitions will only be larger than 
those shown. The dashed lines depict the resulting wall velocities when including the gauge 
boson cubic terms and friction. The values of v w are smaller in this case, although for Set 
2 the gauge boson contribution makes a less significant difference. This suggests that our 
gauge-invariant treatment should provide a reasonable, though rough, estimate of the wall 
velocity when v w is not too large. 

Fig. 2 confirms our intuition from Sec. 6.1: larger thermal masses for the singlet and 
SM-like Higgs field result in slower bubble walls. Larger thermal masses trap the fields 
in the high-temperature minimum and delay the phase transition to lower temperatures. 
This yields smaller changes in the VEVs for a given phase transition strength, and hence a 
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< t>h{T n )/T n 

Figure 2. Wall velocities for the xSM parameter space described in the text. The solid (dashed) 
curves depict the results neglecting (including) the SU(2)l gauge boson contributions to the finite 
temperature effective potential and friction. No subsonic solutions are found with <ph(T n )/T n > 1 
(> 1.1) for the points in Set 1 neglecting (including) the gauge bosons. The curves corresponding 
to Set 2 would extend beyond 4>h(Tn)/T n = 1.1, however the perturbative fluid approximation 
begins to break down significantly for stronger transitions, and so we restrict our results to the 
region shown. The red dotted line shows the speed of sound in the plasma, above which non¬ 
local electroweak baryogenesis is not possible. Note that we have searched exclusively for subsonic 
solutions to the equations of motion. 


smaller pressure difference between the phases. The friction on the bubble wall also tends 
to be enhanced for larger thermal masses. 

Interestingly, for strong first-order phase transitions, we find that subsonic solutions 
to the equations of motion may not exist. This is because as v w —> c s , the background 
temperature contribution begins to dominate in the Higgs and singlet field EOMs (it is 
proportional to l/(Cg — v^)). As pointed out in Ref. [54], the background terms typically 
enter with a relative sign to those from the heavy species, thus reducing the total friction 
for subsonic deflagrations. This behavior is seen for Set 1 in Fig. 2: no subsonic solution 
exists for the gauge-invariant case with <ph{T n )/T n > 1. Including the gauge-dependent 
terms, subsonic solutions can extend up to 4>h(T n )/T n ~ 1.1, but not higher. We conclude 
that viable non-local electroweak baryogenesis in singlet-driven models is incompatible with 
very strong first-order phase transitions, at least in some cases. This can be at odds with 
sphaleron suppression inside the bubble, as seen for Set 1. 

Even if a subsonic solution exists, the bubbles tend to expand rather quickly from 
the standpoint of successful EWB. For example, previous studies of CP-violating sources 
in the MSSM [27-29] suggest that electroweak baryogenesis tends to be most efficient for 
v w ~ 0.01, while Fig. 2 indicates that v w > 0.2 for most points featuring a strongly first- 
order phase transition. Viable bayogenesis in singlet-driven scenarios may thus require 
substantially more CP-violation than in models with slow walls (such as the MSSM with 
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M T n)/Tn 4>h(T n )/T n 


Figure 3. Late-time bubble wall profiles relevant for electroweak baryogenesis obtained by solving 
the wall equations of motion. The solid (dashed) curves depict the results neglecting (including) the 
SU{2)l gauge boson contributions to the finite temperature effective potential and friction. The 
top panel shows the singlet field offset, while the bottom two show the SM-like Higgs and singlet 
wall widths. Bubbles with strong first-order phase transitions tend to feature ~ 5 /T and the 
singlet lagging slightly behind the Higgs field. 


light stops) to overcome the suppression arising from large v w . 

Our methods also allow us to determine the wall widths and offset for the subsonic 
configurations. These quantities are important inputs for microphysical calculations of 
the baryon asymmetry. The resulting bubble wall profiles for Sets 1 and 2 are shown in 
Fig. 3. The offset can change sign, with the singlet held lagging behind that of the SM-like 
Higgs for stronger phase transitions. For <Ph{T n )/T n > 1, the wall widths are typically 
~ 0(h/T). This is substantially smaller than typical values in Standard Model-like cases 
and consistent with the findings of Ref. [48] in the NMSSM. Thin walls follow from the 
large pressure difference due to the changing singlet VEV during the transition. This is in 
fact promising for electroweak baryogenesis, since in many cases the CP-violating sources 
scale as ~ 1 /L w [28, 135]. 

One may ask to what extent we should expect similar results beyond the minimal real 
singlet extension of the Standard Model. After all, the xSM is known to be incomplete 
from the standpoint of electroweak baryogenesis, since it contains no new source of CP- 
violation. However, the model can be modified slightly to incorporate CP-violation by e.g. 
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complexifying the singlet and adding CP-violating Higgs-singlet couplings, or by including 
additional higher dimension CP-violating operators, as in Ref. [112], in cases where the 
singlet VEV vanishes at T = 0. Neither possibility should significantly alter the friction on 
the bubble wall. We expect similar conclusions in other CP-violatiug extensions of the xSM. 
Our findings followed primarily from the form of the friction, which is dominated by the top 
quarks and gauge bosons for the SM-like Higgs field and the singlet and Higgs excitations for 
the singlet field. As long as this is the case, the results beyond the minimal model should be 
qualitatively similar to those we have found here. In fact, this is not unreasonable: it would 
be difficult for new states to couple as strongly to the Higgs field as the top quark without 
violating existing phenomenological constraints, for example. Regardless, the methods and 
ingredients presented in sections 4-5 can be used to determine the wall velocity beyond the 
minimal xSM, although this may require computing additional interaction rates involving 
the excitations of the new species in the plasma. 

7 Summary and Conclusions 

In this study, we have seen how to compute the electroweak bubble wall velocity at singlet- 
driven first-order phase transitions. This extends previous work which applied to the Stan¬ 
dard Model- and MSSM-like cases. For concreteness, we framed our discussion in the real 
singlet extension of the Standard Model, or xSM, although our methods can be used in 
other models involving singlets at the electroweak phase transition. 

Some of the key findings of this study are: 

• As anticipated, bubbles tend to expand rather quickly at first-order phase transitions 
driven by tree-level cubic terms in which the singlet vacuum expectation value changes 
appreciably. We have found v w > 0.2 for all points with (j>h(T n )/T n > 1 in the xSM. 
These wall velocities may be compatible with electroweak baryogenesis in some cases, 
provided a sufficiently strong source of CP-violation. One should bear in mind, 
however, that the free-particle approximation made for the singlet excitations may 
overestimate the corresponding friction, and thus lead to an underestimate of v w . 

• The most promising parameter space for slower bubble walls, and hence for elec¬ 
troweak baryogenesis, features larger thermal masses for the singlet and SM-like Higgs 
field. This translates into larger values of 02 and 64 in the xSM. 

• Strong phase transitions may exhibit no subsonic solution and hence not allow for vi¬ 
able transport-driven electroweak baryogenesis. For example, considering a particular 
set of parameters in the xSM, we found no points with v w < c s for (f>h(T n )/T n > 1.1. 

• A gauge-invariant estimation of the bubble wall velocity is possible and should ap¬ 
proximate the full solution rather well for the slowest walls. The (gauge-dependent) 
SU(2)l gauge boson contributions become important for faster moving bubble walls. 

• Wall widths are typically of order ~ 5/T for strong first-order phase transitions as 
required for electroweak baryogenesis. These values are considerably smaller than 
their Standard Model analogs which are often used in the literature. 
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Along with the above points, our treatment of the friction on the bubble wall can be 
useful in various applications related to the electroweak plasma. For example, the interac¬ 
tion rates computed for the top quarks and Higgs bosons can be used to extract diffusion 
constants for these species, valid at leading log order and including the effects of hard ther¬ 
mal loops, which are important in transport calculations for electroweak baryogenesis. A 
rough estimate along the lines of Ref. [21, 77] suggests D ~ 13/T, Dt ~ 2/T. 

Our results are promising from the standpoint of observable gravitational radiation. 
The peak amplitude for the stochastic gravity wave background produced at a phase tran¬ 
sition is enhanced for faster moving bubbles and larger pressure differences. The ingredients 
presented in this study can be used to more precisely compute the resulting gravity wave 
spectrum in concrete models involving singlets. Recent work [67, 68] suggests that the peak 
amplitude of the signal from a strong electroweak-scale phase transition may in fact be sig¬ 
nificantly larger than previously realized, ft would be interesting to analyze singlet driven 
phase transitions given these new hydrodynamic insights along with our predictions for the 
wall velocity in concrete models. An observable gravity wave signal could provide exciting 
(indirect) evidence for a first-order phase transition, and possibly electroweak baryogenesis, 
in the early Universe. 

There are several ways to improve over the methods presented this study. One might 
hope to move beyond the simple fluid approximation to be able to study stronger phase 
transitions. Also beneficial would be a full leading-order determination of the quasiparticle 
interaction rates entering the Boltzmann equations for the various perturbations. Other 
improvements include accounting for the effects of the spherical bubble geometry on the 
hydrodynamics, formulating a gauge-independent treatment incorporating the gauge and 
Goldstone bosons (which is a difficult problem), and considering full numerical solutions to 
the bubble wall equations of motion rather than utilizing an ansatz. These improvements, 
required for more precise determinations of the wall velocity, would become much more im¬ 
portant if the LHC or a future collider were to unearth direct evidence for a singlet-extended 
Higgs sector. In the interim, we expect our methods to provide a decent approximation of 
the bubble wall velocity in singlet-driven scenarios, which remain a particularly compelling 
setting for electroweak baryogenesis in light of current experimental constraints. 
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A Comparing Effective Interaction Rates 

The effective interaction rates we have computed and listed in Equations 4.25-4.29 differ 
by 0(1) factors from those appearing in the classic references [20, 77]. In this appendix, we 
show that this discrepancy can be explained by the slightly different set of approximations 
made in evaluating the integrals analytically in previous works. Technically, both the 
methods used in Refs. [20, 77] and in this work are valid ‘leading log approximations’ to 
the full first-order result, in that only processes contributing logarithmically to the collision 
integrals are considered. They differ primarily in their treatments of the non-logarithmic 
pieces of the various momentum integrals. The discrepancies can therefore be understood 
to demonstrate the uncertainties associated with the leading log approximation and the 
importance of performing a full leading-order calculation for more precise results in future 
studies. 

To understand the different approaches to evaluating the collision integrals, let us con¬ 
sider the process tt —> gg. First of all, Ref. [77] does not include the symmetry factor for 
the corresponding matrix elements, resulting in a factor of two discrepancy before evaluat¬ 
ing any integrals [121], Including the symmetry factor, the leading-log matrix element is 
~ 64/9g r |u/t, where we have averaged over the top quark degrees of freedom. The result¬ 
ing contribution to r* (lt , evaluating the integral numerically and including the top quark 
momentum-dependent self-energy, is 

Ar^ M « 1.1 X io- 3 t, (at) 


whereas Ref. [77] reports 

A dM = ^xflo 6 ^T»3.8xlO-r. (A.2) 

including the correct symmetry factor. Starting from the same matrix element, the different 
methods for evaluating the integrals results in almost a factor of 4 difference in the result. 

Simply evaluating the integral in Ref. [77] numerically, but including the thermal mass 
in the propagator instead of the HTL momentum-dependent self-energy, we find 

Ar£ 1)t « 1.5 x 10- 3 T, (A.3) 


suggesting that the simple propagator replacement over-estimates the integral by about 
40%, but does not account for the whole discrepancy. In fact, most of the difference comes 
from the various approximations made to arrive at the analytic result in Eq. A.2. In 
particular, all non-logarithmic contributions are dropped, while numerically evaluating the 
integrals includes all of the various contributions. 

For example, the final result for AT^ f in Ref. [77] contains an integral over the plasma 
frame angle 9 between p and k, 


dcosO- log 


2 |p| |k| (1 — cos 9) 


mt 


-1 + log 


4 |p| |k 

9 

mf 


(A.4) 
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Moore and Prokopec drop the constant piece, keeping only the logarithm. However, the 
contribution of the constant piece is numerically comparable, and of opposite sign, to the 
logarithmic term. Performing the remaining integrals over |p|, |k|, the contribution without 
the constant term is ~ 3.1 X 10~ 3 T, while including it yields ~ 1.8 X 10 _3 T, which is 
significantly closer to the results we have obtained. Similar approximations are made in 
performing the other integrals, and for the other rates. 

It is worth reiterating that neither approach includes all processes contributing at 
leading order in the gauge couplings. Our interaction rates should simply be viewed as a 
slightly different approximation to the full leading order result. Future studies of the wall 
velocity in different models may find it beneficial to compare the results obtained from our 
interaction rates and those of Ref. [77] to assess the uncertainty expected from the leading 
log approximation. 
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